Case1 - University of Nevada

108024502 林書宇 , 108024507 張文騰 , 108024701 張珮榕 , 108024702 陳信諺

Abstract

此資料是由University of Nevada的梁靜農教授所提供。資料中涵蓋17個人(14位正常人及3位中風患者)其脛前肌與比目魚肌抑制關係的實驗相關結果，而資料分析目的在於想要了解正常人與中風患者在脛前肌與比目魚肌抑制關係的差異。在分析結果中，我們使用一個二階模型來描述資料的34隻腳，且用 Euclidean Distance 針對正常人/中風患者進行分類。最後，我們嘗試針對抑制比提出一個新的想法，並且與目前公認的抑制比公式進行比較。

Problem Formulation

Problem 1

找出比目魚肌的出力 \(Y\) 與脛前肌的出力 \(X\) 間的關係

針對資料中的34隻腳，我們目標在於建立一個模型，讓每個人的每隻腳都能使用該模型去描述“脛前肌收縮狀態”對“比目魚肌收縮狀態”的影響，且模型能夠對於每隻腳都有不錯、合理的描述。

若以統計模型表示，則此問題可以描述成: \(比目魚肌收縮程度 = f(脛前肌收縮程度)+\epsilon\)

  

Problem 2

找出可以區分正常人、中風人(好、壞腳)的關鍵特徵

自Problem 1的模型中，希望可以從資料找到一些特徵，且該特徵在中風患者的表現與正常人的表現上，差異能夠越大越好。

而此處的特徵是指模型的參數估計量、以及評估參數估計量好壞的數值(如:標準誤、P-value)等。

若以數學的概念表示，則此問題可以描述成:
尋找\(f(X)\)的特徵, \(w_{i,j} \ni w_{t_1,j} \in W_1 , w_{t_2,j} \in W_2\) ,   

where

\((1) i = 1 , ... , 34 , j = 第j種特徵\)
\((2) W_1: 正常人的特徵 , W_2:中風患者的特徵 ( i_1 \in {1,...,28} , i_2 \in {29,...,34} )\)
目標 : \(w_{t_1,j}與w_{t_2,j}的差異越大越好\)

  

Problem 3

將正常人、中風人好腳、中風人壞腳做分類

從Problem 2所找出的特徵中，使用可以清楚區分中風患者與正常人的特徵制定出一個標準，以便將每隻腳對於是否屬於中風患者而進行分類。

(若結果允許，此處也可能考慮將資料中的34隻腳針對這些特徵分成三類，分別為: 中風患者受中風影響的腳、中風患者的正常腳、正常人的腳)

  

Problem 4

抑制比的估計法

對於目前公認的抑制比估計方法，將Problem 1 的模型進行套用與修改。

其概念為 : \(\frac{f(X_{max}) - f(X_{min})}{f(X_{max})}\)，其中f(.)為Problem 1所使用的函數，\(X_{max}\)以及\(X_{min}\)則為每隻腳的最大/最小\(X\)值。

  

或是嘗試提出一個新的估計量，並且與原本的抑制比估計方法進行比較。

其概念為 : \(脛前肌對比目魚肌的抑制比 = g \left(f(.),X_{max},X_{min}\right)\)，其中\(f(.)、X_{max}\)以及\(X_{min}\)如上方所述。

  

Data Cleaning and Processing

資料處理示意圖如下:

其中:

  1. 人編號: 1-14為正常人，15-17為中風患者。編排順序則與原先所得資料相同 (NI001, …, NI016, CVA003, CVA004, CVA005)。
      後續我們會以“第幾個人”表示之。其中第15至第17人分別為中度/重度/輕度中風患者。

  2. 左右/好壞腳: 在正常人的資料中，可以分為左右腳。而中風患者並無左右腳的資訊，但有好壞腳的資訊。其中好腳為中風患者的正常腳，壞腳則為中風患者受中風所影響的腳。

  3. 腳編號: 資料中34隻腳的編排順序與人編號相同。而正常人以先左腳後右腳的順序編排，中風患者則以先好腳後壞腳的順序編排。
      (NI001左腳, NI001右腳, …, NI016左腳, NI016右腳, CVA003好腳, CVA003壞腳, …, CVA005好腳, CVA005壞腳)

  4. scale_X 的數值為: AntagonistBG (X)的數值，針對每一隻腳，分別除以該腳的最大AntagonistBG (X)數值。

  5. 左右腳、好壞腳、沒病/中風、ID2-14,16, 17: 簡述之，由於我們是使用一個模型涵蓋這34隻腳，以進行模型配適及描述狀況。這些變數不妨可以看成是可以從模型中辨別是哪隻腳的變數。若以統計的角度進行說明，由於資料是nesting的結構，因此我們須先將資料區分成正常人以及中風患者，其中正常人又可以分成左腳與右腳、中風患者又可以分成好腳與壞腳。而關於人的區分方式，則是以Dummy Variable(ID2-14,16, 17)處理。

Exploratory Data Analysis

畫出1-14個正常的左右腳圖形，以下這個兩個圖的X軸我們放的是已經Scale過的X，Y軸放的是ppresponse，從左上角的第一張分別是第一個正常人的左腳和右腳，再來是放第二個人的左右腳以此類推，第二張圖的最後三個圖為三個病人分別為中度重度以及輕度，那紅色點代表病人的好腳，綠色點代表病人的壞腳。

可以發現在正常人中每個人的圖形有一定的規律，大致為一個圓滑的L型圖案，並且每個人的轉折點位置不一致，有的人的轉折點比較靠近Scale X偏小的地方，有的比較靠近Scale X，並且有些人的左右腳其實仍有差距存在，但在病人的雙腳中，可以發現病人壞腳較無Ｌ的規律，但病人好腳中仍有存在著Ｌ的規律，但較不明顯。

再來我們好奇，假如把同一個人的資料只分成休息與沒休息時圖形呈現為何，以下兩張圖的X軸為ScaleX，y軸為ppResponse，並且每個人的左右腳疊加在一起，紅色點為休息狀態的資料，綠色點為出力狀態的資料點，可以發現不管左腳右腳好腳壞腳的休息狀態都會在前端，因此可以聯想到抑制比，每個人的休息狀態的點，可以當成每個人的Baseline，因為那是每個人的最原始的狀態，再和出力狀態的做比較和轉換，可以換得每個人的中風狀態。

假如我們把所有人的所有數據畫在同一張圖上，x軸放ScaleX，y軸放ppresponse，紅色的點代表所有正常人的數據點，綠色點代表所有病人的數據點，由這張圖可以確認，病人的規律明顯和正常人不同。

緊接著我們倆探討每一隻腿之間的差距關係，因此我們畫了統計上常用的Box plot。 Box plot的中文翻譯是箱型圖，他是一種可以探討數據分散情況的統計圖，箱子的中間那條線代表的是全部數據最中間的那個數，也就是前50%的資料點，統計上稱為「中位數」(median)，而箱子的下方那條線代表的是前25%的數據量，那箱子的上方那條線代表前75%的數據量，箱子最上以及最下方延伸出去的線則可以代表說這筆數據有多少離群點，因此我們可以用箱型圖來簡單探討數據的分佈。 下面這張圖是全部人的腳，總共34隻，可以觀察到就算是同一個人，他們的兩隻腳也可能有差距，像是第一位正常人的兩隻腳，中位數的地方就有明顯的差距，重度患者的雙腳差距甚大，幾乎所有數據的分佈位置都不一樣，也就是說從這張圖，我們覺得每隻腳的差異大；此外，我們也觀察出來左右兩隻腳、好腳與壞腳的箱型圖也有些差距，因此我們偏向對每一隻腳做建模的動作，而不是對每隻個人。

因為從上一張整體圖形，我們發現可能病人中的好腳與壞腳之間有差距，所以我們對所有病人的好腳壞腳綜合起來看，紅色代表是病人的壞腳，綠色代表病人的好腳，可以發現病人的壞腳中位數都比好腳的中位數低，第一四分位距跟第三四分位距都相差有一段距離。

接著看正常人中的左右腳是否有差異，我們把所有的正常人左右腳的數據綜合來看，紅色箱型圖為正常人的左腳，綠色為正常人的右腳，y軸為ppresponse，由此圖似乎沒有像病人的好壞腳差距那麼大，在正常人的左右腳中，中位數都落在差不多的位置，第一以及第三四分位距也相近，唯一左腳的離群點比較多些，推論應該是屬於第一位正常的腳，因為第一位病人的腳ppresponse有偏大的值。

接下來我們想用一條平滑曲線來代表一隻腿的圖形，那這條平滑曲線我們會採用統計母數中的一個方法稱為「局部多項式」，首先我們會先開一個小的區間，在這個區間中用數學泰勒展開的方式去逼近這個區間中的函數，所以下兩張圖是我們用此方法去畫出每一隻腳的一條平滑函數，接著的內容會去探討這些函數彼此間的關係。

首先我們先把所有正常的左右腳曲線分別用黑色與紅色來代表他，綠色線代表病人的好腳加入，下圖的x軸為Scale X，y軸放ppresponse，我們可以明顯地看出這三個顏色的曲線的圖形是彎曲曲線，越靠近x偏小的地方正常人的腳彎曲程度比較大，但病人的好腳的彎曲程度明顯比正常人來的低。

假如單看加入病人壞腳的圖形進來看的話，可以發現病人壞腳的曲線幾乎沒有像是正常人好腳一樣的彎折狀況，比較偏向水平線。

那接下來我們對這34條線取一次微分，在數學上一次微分的意思代表找線的每一個點的斜率，那我們一開始想找斜率的想法來自說正常人的每一條線都有一定的曲度，可是病人的雙腳的曲線的曲度卻沒有正常人來的彎曲，那曲度可以用斜率來表達，因此我們想看看到底差距有多少。 下面的圖的x軸為Scale x，y軸為估計ppresponse的斜率值，我們先畫出正常人的28隻腳的曲線的一次微分曲線，可以發現在x偏小的地方大家的斜率都由負的轉正的，也就代表說一開始大家會遇到斜率很陡的開頭，然後漸漸地轉彎後變平緩。

假如我們把病人的雙腳曲線也加入探討的話，綠色曲線為病人的好腳，藍色曲線為病人的壞腳，可以很明顯的發現病人的斜率真的如同我們想的一樣十分的平緩，與正常人的雙腳曲線有明顯的差異存在。

Model Building

二次多項式回歸模型

在這部分，我們推薦的模型為二次多項式回歸模型，其具如下形式：

其中有幾個地方須特別注意：

  1. 每隻腳的Y皆做了 \(\log\) 轉換
  2. 每隻腳的X皆做了 \(\frac{\text{antagonist X}}{\max (\text{antagonist X})}轉換\)，分母的\(\max\)為原始資料中特別提供的最大值
  3. 殘差並非為一般所見 \(\epsilon\overset{i.i.d}{\sim} N(0,\sigma^2 I)\) 同質變異數之形式，為了使模型配適的更好，我們多考慮了權重，而權重反映在殘差上後即有此種假設(異質變異數)，\(\epsilon\overset{i.i.d}{\sim} N(0,\sigma^2\Sigma)\)，其中 \(\Sigma = \begin{bmatrix} \frac{1}{x_1^2} & \cdots & 0\\ \vdots & \ddots \\ 0& & \frac{1}{x_n^2} \end{bmatrix}\)
  4. 模型中粗體係數含有15個係數
  5. 此模型總共有 \(102\) 個變量(含截距)，且一共有 \(4462\) 個樣本

模型的推演過程

  1. 為什麼不單純使用 \(E[\log(Y)] = \beta_0 + \beta_1X +\beta_2 X^2\)，而加入很多變量使模型看起來如此複雜？
  2. 為什麼反映變數是 \(\log(Y)\) 而非 \(Y\)？
  3. 為什麼要考慮權重 (也就是殘差有異質變異數的假設)？

上述所列 \(1,2,3\) 分別對應下圖中相應的部分

1. 為什麼不單純使用 \(E[\log(Y)] = \beta_0 + \beta_1X +\beta_2 X^2\) 去配適每一隻腳？

下圖中 \(A\) 為我們所推薦輸入為整筆資料 \(4462\) 個樣本的方法，\(B\) 為每隻腳分別去配適的方法，輸入為每隻腳分別的樣本數

  1. \(A\) 方法可以得到統一的判定係數 \(R^2\)
    • 分開每隻腳去配適(即用 \(B\) 方法)，我們會得到 \(2 \times 17\) 個 \(R^2\)，用 \(A\) 方法可以得到一個所有腳的overall \(R^2\)；如此我們可以更加方便的比較各個模型的好壞
  2. \(A\) 方法中p-value的基準相同，\(B\) 方法中不同隻腳所得的係數估計之顯著程度，不能跨腳比較
    • 分開每隻腳去配適(即用 \(B\) 方法)，我們會得到 \(2 \times 17\) 個不同的 \(\hat{\sigma}^2\) (因為每隻腳的 \(\hat{\sigma}^2\) 會分別用其殘差運算所得的 \(MSE\) 去作估計)，而此 \(MSE\) 的大小將直接影響信賴區間的寬度，進而影響p-value的顯著程度，若我們採取 \(A\) 方法，我們會得到唯一的 \(\hat{\sigma}^2\)，模型中所有的區間是統一的 \(\hat{\sigma}^2\) 運算而得，故係數間可互相比較
  3. 兩個模型可以算出相同結果，但 \(A\) 方法係數可以讓我們了解各個變數的效應
    • 將於以下介紹

接下來我們將以 \(2\) 號正常人的左腳來做舉例：

  1. 上圖中於最上方資料集的部分，橘色底即為 \(2\) 號正常人的左腳之變數輸入值 (\(2\) 號左腳總共有\(114\)筆資料，資料間除了 \(\text{scale_X}\), \(\text{ppResp} (Y)\) 不同外，其餘column都相同如下所示 (\(ID4\) 後的column因版面關係故省略，可參考上方橘色底的部分除了 \(ID2\) 為 \(-1\) 外其餘皆為 \(0\)))
Table 1: 2號正常人的左腳
人編號 左腳-1右腳1 好腳-1壞腳1 ppResp (Y) 沒病0中風1 scale_X ID2 ID3 ID4
2 -1 0 2.941895 0 0.2157098 -1 0 0
2 -1 0 2.973498 0 0.2232633 -1 0 0
2 -1 0 2.928929 0 0.2237095 -1 0 0
2 -1 0 2.830877 0 0.2255610 -1 0 0
2 -1 0 2.944812 0 0.2293764 -1 0 0
2 -1 0 2.885332 0 0.2294466 -1 0 0
  1. 了解了輸入值後，我們可以看到資料集下方 \(A\) 模型藍色底的輸出部分，我們將挑出一些係數來做解釋：
    • \(\beta_{左腳-1右腳1} = -1.1188911\)： 這項係數代表只要是左腳，平均而言 \(\log(Y)\) 都會增加 \((-1) \times (-1.1188911)\)
    • \(\beta_{scale\text{_}X:ID2:左腳-1右腳1}\;\;\;= 0.2835283\)： 這項係數代表只要是左腳且 \(ID = 2\)，平均而言 \(\log(Y)\) 會增加 \((-1) \times (-1) \times 0.2835283 \times X\)
    • \(\beta_{沒病0中風1:I(scale\text{_}X\text{^}2)} = -10.7744410\)： 這項係數代表只要是沒病，\(\log(Y)\) 不會改變，因為 \(0 \times 10.7744410 \times X^2 = 0\)
  2. 接下來我們可以看到上圖中模型輸出的部分，我們發現到不論是截距項係數、一次項係數、二次項係數，\(A\) 方法和 \(B\) 方法皆有相同結果，而 \(A\) 方法成功將 \(B\) 方法中的係數拆解成多個效應，我們因此可以知道說到底是哪些效應(左右腳、好壞腳、有無中風)的影響促成截距項係數、一次項係數、二次項係數的改變。

2. 為什麼反應變數是 \(\log(Y)\) 而非 \(Y\)？

  1. 因為用 \(Y\) 作為反應變數會違反模型假設 (\(Y\)和殘差不為常態分配)，而 \(\log(Y)\) 則較符合假設
  2. 因為 \(\log(Y)\) 配適出來的圖形較合理，原因如下：
    • 尾端向上較輕 (比目魚肌出的力應與脛前肌出的力成反比)
    • \(Y\) 值不會發生小於\(0\)的狀況 (比目魚肌出的力不應該低於\(0\))

細節如下：

  1. 關於上方理由第一點，我們可以由下方圖發現 \(Y\) 值呈現嚴重右偏分配，故須做 \((\log)\) 轉換

  1. 關於第二點，首先我們可以先觀察如果只用 \(Y\) 而非 \(\log (Y)\) 配適出來的結果，我們可以發現到尾端嚴重向上、且比目魚肌出力有低於\(0\)的狀況發生

  1. 接著我們可以觀察如果使用 \(\log (Y)\) 配適出來的結果，我們可以發現到尾端向上趨勢稍微緩解、且比目魚肌出力已不再有低於\(0\)的狀況發生

3. 為什麼要考慮權重？

考慮權重之所以可以減輕尾端向上情況的解釋如下：

  1. 首先，我們考慮一個只有 \(5\) 個觀測值的資料，圖形如下所示，紅色線為納入權重前，而藍色線為納入權重後；我們可以發現到藍色線較紅色線靠近最後一點，原因如下

  1. 原始 \(RSS\text{(Residual sum of Squares)}\)有如下形式： \[RSS = \sum_{i=1}^n(Y_i-\hat{Y_i})^2 = \sum_{i=1}^n(Y_i-x_i\hat{\beta})^2\] 考慮權重後的 \(RSS\) 有如下形式 (令權重 \(w_i\) 為 \(x_i\) )： \[RSS = \sum_{i=1}^nw_i(Y_i-\hat{Y_i})^2 = \sum_{i=1}^nw_i(Y_i-x_i\hat{\beta})^2 = \sum_{i=1}^nx_i(Y_i-x_i\hat{\beta})^2\] 我們可以發現到，當 \(x_i\) 越大時，\(x_i(Y_i-x_i\hat{\beta})^2\) 就會越大，而我們在找迴歸線時是透過最小化\(RSS\)的方式去尋找的，也就是說，越大\(x_i\)值的地方我們要配適越好，不然會被懲罰(權重 \(x_i\))太重。故考慮權重 \(w_i = x_i\) 的迴歸線會對越大的 \(x_i\) 配適越好。因此可以成功將尾端壓下來。

結果如下，我們已成功將尾端壓下來且前端配適的仍可接受，此時的 \(R^2 = 0.8564\)

Analysis Result

Problem 1

根據以上分析，我們認為 \(Y\) 與 \(X\) 直接的關係比較難去model，故我們透過將 \(Y\) 做 \(\log\) 轉換來描述他們之間的關係。

Problem 2

我們可以發現到壞腳的曲線都很不明顯，且對於每隻腳斜率看起來也不盡相同，故我們考慮將截距項係數、一次項係數及二次項係數拿出來看看，是否在各種腳中看得出不同的趨勢

Problem 3

Euclidean Distance 的公式為： \[D = \sqrt{r_1^2 + r_2^2 + r_3^2}\] 其中 \(r_1, r_2, r_3\)分別代表截距項係數、一次項係數及二次項係數

結果 \((34 \text{隻腳的}D)\) 呈現如下圖：

  1. 此方法有成功將中風人(紅、橘點)與正常人(藍、深藍)點分得很開
  2. 對於壞人好腳(橘點)與壞人壞腳(紅點)的分類，我們並沒有成功將其完全分開，可以發現到壞人壞腳(紅點)的值應該要偏低(透過觀察最重症病人 \(16\) 的紅橘點分布)，但是第 \(17\) 號病人他的紅點卻比橘點略高

Problem4

由於原先對於抑制比的定義偏重於將脛前肌於放鬆狀態的比目魚肌與脛前肌在最大收縮狀態時的比目魚肌進行比較，在這樣的定義下由於只注重左、右兩端點，造成固定兩端點下，其實不同的曲線都會得到同樣的抑制比。對於在脛前肌較放鬆時，比目魚肌如果很快地被迫放鬆(曲線很陡)，那這樣的抑制效果應該是比慢慢地被迫放鬆(曲線較平緩)來的強大。因此我們將嘗試於此節提出新估計量，試圖將兩端點之外的資訊也包含在其中(包括斜率、彎曲程度等等)，並於後續討論其優缺點及目前的缺陷，再與原本定義的抑制比進行比較。

對於每一隻腳，我們可以根據前面線性模型的分析結果繪製一條曲線代表之，我們若將此曲線的最左端與最右端連線將可形成一條直線，再將此腳的脛前肌最放鬆X值、最收縮X值分別繪製一條垂直線，我們便可將整個圖形區分為A、B兩個部分，如下圖所示：

我們可將新的估計量定義為

\[\frac{A面積}{A面積+B面積}\]

假設上方這條曲線往右上角靠近，這就造成了A面積縮小，B面積變大，A面積+B面積不變，但估計量的估計值下降。變動後的曲線直觀上的意義為抑制作用的效果降低了，沒有像原先那麼強大，這就會與估計值的表現一致；若曲線向左下角靠近，A面積變大，B面積變小，估計量的估計值增加。變動後的曲線直觀上的意義為抑制作用的效果增強了，脛前肌在收縮程度偏小的地方只要增加一點點收縮力，比目魚肌就會受到抑制作用非常大的影響而快速放鬆，這依舊與估計值的表現一致。

假設曲線向上平移，也就是脛前肌出力時，比目魚肌的放鬆狀況不理想，抑制作用在脛前肌較用力時效果不顯著，則A面積不變，B面積增加，估計量的估計值下降，值的表現與解釋一致，反之亦然。

依據上面的定義方式，新估計量的估計結果如上，中風患者其兩腳的數值值皆小於0.46。

雖然上面這個新估計量的定義方式看起來還不錯，有將曲線的斜率以及彎曲程度表現於估計量中，但我們會發現對於一條被估計為“直線”的腳，其在此定義下的估計值會變成0，這並不是我們所樂見的，因為即使該腳的曲線被我們估計為直線，其抑制效果依舊正常作用，被我們估計為0是非常不合理的，因此我們可將這個新估計量的定義方式做些修改：

我們一樣延續原本的定義概念，但將“A面積的定義”改動一下成為

在這樣的定義下，原先提到的“斜直線估計”問題，就可以被我們解決了，被我們估計成斜直線的腳，其新定義估計量的估計值會成為0.5。

我們若使用修正過的定義法去估計34隻腳，其結果如下：

可以發現第15位(中度中風患者)的好腳其估計值離正常人很接近，而第17位(輕度中風患者)的好腳估計數值竟比其壞腳低了一些。雖然此修正估計法解決了斜直線的矛盾問題，但對於每隻腳的估計效果似乎沒有修正前那麼好。

接下來，我們將列出3種數值估計法的表格進行比較，第一種為我們的新定義估計法(未修正)，第二種為我們的新定義估計法(修正後)，第三種為梁老師原始的估計法

腳編號 新定義估計法.未修正. 新定義估計法.修正後. 原抑制比 未修正的新定義與原抑制比的差 修正後的新定義與原抑制比的差
第一人左腳 0.6643591 0.8205890 0.9593923 0.2950333 0.1388034
右腳 0.6853330 0.8211426 0.9202813 0.2349483 0.0991387
第二人左腳 0.6001345 0.7959202 0.9476877 0.3475532 0.1517675
右腳 0.3646117 0.6380014 0.8251990 0.4605873 0.1871976
第三人左腳 0.8080988 0.8922559 0.9852628 0.1771639 0.0930069
右腳 0.7723058 0.8796734 0.9541657 0.1818599 0.0744923
第四人左腳 0.7645952 0.8448807 0.9637909 0.1991957 0.1189102
右腳 0.7507638 0.8694631 0.9528360 0.2020722 0.0833729
第五人左腳 0.7699699 0.8797585 0.9687903 0.1988204 0.0890318
右腳 0.6697784 0.8181601 0.9451243 0.2753459 0.1269642
第六人左腳 0.6054163 0.7620144 0.7922055 0.1867892 0.0301911
右腳 0.6160788 0.7250051 0.8433736 0.2272948 0.1183685
第七人左腳 0.6036473 0.7820754 0.8159016 0.2122544 0.0338262
右腳 0.6045930 0.7812861 0.7574243 0.1528314 -0.0238618
第八人左腳 0.7696326 0.8736498 0.9427268 0.1730942 0.0690771
右腳 0.8148866 0.8981697 0.9472274 0.1323408 0.0490577
第九人左腳 0.7463432 0.8580276 0.9069761 0.1606329 0.0489485
右腳 0.7594782 0.8665561 0.9425241 0.1830460 0.0759680
第十人左腳 0.4703902 0.5918360 0.6408847 0.1704945 0.0490487
右腳 0.6361387 0.7827567 0.9195205 0.2833818 0.1367638
第十一人左腳 0.7210669 0.8437953 0.9435827 0.2225158 0.0997873
右腳 0.6608476 0.7907395 0.9337314 0.2728839 0.1429919
第十二人左腳 0.7926802 0.8928860 0.9127581 0.1200779 0.0198721
右腳 0.7096743 0.8488205 0.9220002 0.2123259 0.0731798
第十三人左腳 0.5621696 0.7351018 0.7100528 0.1478832 -0.0250490
右腳 0.6694438 0.8177545 0.9707575 0.3013137 0.1530030
第十四人左腳 0.6892286 0.8306041 0.9173961 0.2281676 0.0867921
右腳 0.6073787 0.7880062 0.9581669 0.3507883 0.1701608
第十五人好腳 0.4562632 0.6849085 0.8160320 0.3597687 0.1311235
壞腳 0.2044259 0.4676313 0.5017589 0.2973330 0.0341276
第十六人好腳 0.3970001 0.5546819 0.7550105 0.3580104 0.2003285
壞腳 -0.0031033 0.0435435 0.3560340 0.3591373 0.3124905
第十七人好腳 0.4223822 0.5701120 0.8139385 0.3915563 0.2438265
壞腳 0.4580043 0.6277831 0.7638456 0.3058413 0.1360625

在比較新估計結果與原估計結果時，新定義(未修正)與原定義差距超過0.27的腳且新定義(修正後)與原定義差距超過0.12的腳有：第1人的左腳、第2人的雙腳、第5人的右腳、第10人的右腳、第11人的右腳、第13人的右腳、第14人的右腳以及15、16、17人的雙腳。 可以發現原先的抑制比估計(橘色)對於正常人與中風好腳來說，估計值十分相近，特別偏低的估計僅有中度中風患者與重度中風患者的壞腳。但以我們定義的未修正新估計(黃色)方式，除了中風患者的壞腳外，連好腳的估計值也可以被我們壓下來了，但第二人的右腳與第十人的左腳也是有偏低的狀況。使用修正後的新估計方式(綠色)雖然解決了直線矛盾問題，但整體的估計卻沒有比較好的區分效果，可能不能很好的表現抑制的概念。

Information summary

我們可以從一開始的探索性資料分析與接下來的線性模型分析歸納出4個重點：

1.從模型觀察出的正常人與中風患者差異

正常人脛前肌與比目魚肌的抑制作用在脛前肌沒有那麼用力的時候會非常強烈(在X值偏小的時候，隨著X值的增加，Y會下降地非常快速)；在脛前肌越來越用力後，抑制作用造成的效果會稍微趨緩些，爾後飽和。也就是當X值繼續增大後，我們會發現Y值的下降速度沒有一開始那麼快，抑制效果開始有了邊際效應，不斷遞減直到飽和。

中風患者的抑制關係不論好腳或壞腳，在脛前肌沒有那麼用力的時候就比正常人低了許多(在X值偏小的時候，隨著X值的增加，Y值的下降比正常人慢了許多)；脛前肌繼續用力的抑制效果變化量與正常人相比也低了很多，正常來說這時候應該要非常快受到邊際效應的遞減，但中風患者的邊際效應沒有這麼強烈，圖形看起來會“比較沒有那麼彎”。中風患者比目魚肌的平均收縮大小低於正常人許多，甚至其分布範圍本來就比大部分正常人低了許多。

2.足以用來區分病人與正常人的特徵

我們於線性模型的分析中，得到了三個可以明顯區分正常人與病人的特徵：截距項係數、scaleX係數、scaleX平方的係數。這三個係數在統計數學上的解釋可以分別對應到Y值平均大小、估計曲線的切線斜率、估計曲線的轉彎程度，這三個特徵從一開始的探索性資料分析就已經出現了。而他們在人類肌肉上的解釋是比目魚肌的平均收縮大小、抑制效果的強度表現、抑制效果隨著脛前肌繼續收縮的邊際效應。

3.中風的好壞程度表現於中風患者的好腳

我們於探索性資料分析中，曾看到中風患者的好腳其實會根據其中風程度而有些微不同的斜率和曲度，整條曲線越像正常人其中風程度就會越低。在線性模型的分析中，我們也看到了中風的程度會表現於中風患者的好腳參數估計值與正常腳的估計值距離。

4.中風人好壞腳的同化現象

中風患者兩腳的平滑曲線很接近，與正常人的表現皆不同，有點像是中風患者的壞腳將其好腳的表現從正常腳拉至壞腳。在線性模型的參數估計也有類似的表現，好腳的各個估計值都偏離正常腳許多，與其壞腳很接近。

Appendix

1.探索性資料分析中的平滑化方法－局部線性回歸

假設二元隨機樣本\((X_1,Y_1),\ldots,(X_n,Y_n)\)來自以下模型：

\[\begin{align} Y = m(X)+\varepsilon, \end{align}\]

其中，Y為連續反應變數，X為連續解釋變數，\(\varepsilon\)為誤差且滿足\(E(\varepsilon)=0\)與\(Var(\varepsilon)=\sigma^2\)，X與\(\varepsilon\)相互獨立。

假設\(m(X_i)\)可以使用泰勒展開式進行逼近

\[\begin{align} m_p(X_i)=m(x)+m'(x)(X_i-x)+ \dots +m^{(p)}(x)(X_i-x)^p/p!, \end{align}\]

對於與x足夠近的\(X_i\)，我們便可以估計\(m(x)=E(Y|X=x)\)及其導數\(m^{(v)}(x)\)。

所謂的局部多項式就是使用加權最小平方(weighted least squares)來局部擬合p階多項式\(m_p(X_i)=\sum_{k=0}^p \beta_{k}(X_i-x)^k\)藉由最小化：

\[\begin{align} \label{WLSS} \sum_{i=1}^{n}[Y_i-m_p(X_i)]^2K_h(X_i-x) = (\textbf{y}-\textbf{X}_p\beta)^T\textbf{W}_x(\textbf{y}-\textbf{X}_p\beta), \end{align}\]

其中，\(y=(Y_1,\dots,Y_n)^T\)，\(\textbf{X}_p = [\textbf{1} \quad \textbf{x}_{1x}\quad \dots \quad \textbf{x}_{px}]_{n\times(p+1)}\)為設計矩陣且\(\textbf{1} = (1,\dots,1)^T\)長度為n的向量，\(\textbf{x}_{kx}=((X_1-x)^k,\dots,(X_n-x)^k)^T\)，\(k=1\dots,p\)，\(\beta = (\beta_0,\dots,\beta_p)^T\)與z有關，\(\textbf{W}_z=diag\{K_h(X_1-x),K_h(X_2-x),\dots,K_h(X_n-x)\}\)其中的\(K(\cdot)\)為非負且對稱的機率密度函數，h為平滑參數，\(K_h(\cdot)=K(\cdot/h)/h\)。

我們可將滿足加權最小平方的解寫成

\(\widehat{\beta}=(\widehat{\beta_0},\dots,\widehat{\beta_p})^T = (\textbf{X}_p^T\textbf{W}_x \textbf{X}_p)^{-1} \textbf{X}_p^T \textbf{W}_x \textbf{y}\)，利用 \(\widehat{m}(x)=\widehat{\beta_0}\) 來估計 \(m(x)\)，用 \(\widehat{m}_p^{(v)}(x)=v!\widehat{\beta}_v\) 來估計 \(m^{(v)}(x)\) (Fan and Gijbels 1996)。

至此，我們可將p階局部多項式的估計式\(m^{(v)}(x)\)寫成矩陣形式(使用與Aurore Delaigle, Jianqing Fan & Raymond J. Carroll 2009同樣的符號表示)：

\[\begin{align}　\label{mpd} \widehat{m}_p^{(v)}(x) = v!\textbf{e}_{v+1}^T(\textbf{X}^T \textbf{K} \textbf{X})^{-1} \textbf{X}^T \textbf{K} \textbf{y} \end{align}\]

其中，\(\textbf{e}_{v+1}^T = (0,\dots,1,\dots,0)\)表示除了第(v+1)元素為1其餘元素皆為0的向量，\(\textbf{y}^T=(Y_1,\dots,Y_n)\)，\(\textbf{X}=\{ (X_i-x)^j \}_{1\leq i \leq n,0 \leq j \leq p}\)，\(\textbf{K} = diag\{ K_h(X_j-x) \}\)。